library(DESeq2)
library(tidyverse)
library(tximport)

tx2gene <-
  read_delim('~/append-ssd/nextflowing/rnaseq_wu2021_batch1/salmon/tx2gene.tsv',
             col_names = F)

wu21_b1 <-
list.files('~/append-ssd/nextflowing/rnaseq_wu2021_batch1/salmon/',
           'quant.sf', full.names = T, recursive = T)

wu21_b1_txi <- wu21_b1 |>
  str_extract('P\\d+') |>
  set_names(wu21_b1, nm = _) |>
  tximport(type = 'salmon', tx2gene = tx2gene)

head(wu21_b1_txi$counts)

sample_info <- read_tsv('Archive/covid19/data/wu-chen-samples.txt')

coldata_b1 <- sample_info |>
  filter(id1 %in% colnames(wu21_b1_txi$counts))

dds <- DESeqDataSetFromTximport(wu21_b1_txi, colData = coldata_b1)

wu21_b1_txi|> glimpse()

data <- read_csv('Archive/covid19/data/wu-chen-salmon.count.csv.gz')

data <- data |>
  select(gene_name, value, id1) |>
  pivot_wider(names_from = id1, values_from = value, values_fn = sum)

data |>
  mutate(across(-1, round)) |>
  write_csv('Archive/covid19/data/wu-chen-RNA-round-count.csv')

data <- read_csv('Archive/covid19/data/wu-chen-RNA-round-count.csv')

data <- data |>
  column_to_rownames('gene_name')

subset_info <- sample_info %>%
  filter(id1 %in% colnames(data)) |>
  column_to_rownames('id1')

subset_info <- subset_info |>
  mutate(severe_bool = str_starts(clinic, '3|4'))

dds <- data |>
  DESeqDataSetFromMatrix(
  colData = subset_info,
  design = ~ severe_bool)

# pre-filtering data of too low counts
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

# relevel to specify the control group (reference)
levels(dds$severe_bool)
dds$clinic <- relevel(dds$clinic, ref = "1_asymptomatic")

# rlog() may take a long time with 50 or more samples,
# vst() is a much faster transformation

# set parallel workers reduce time cost from 188s to 96s for 63*24200 matrix
dds <- DESeq(dds, parallel = TRUE)
res <- results(dds)
summary(res)

deg_tibble <- res |>
  as_tibble(rownames = 'gene')

deg_tibble |>
  filter(gene == 'TRPM2') |>
  mutate(logFC = log2FoldChange, adj.P.Val = padj,
         dataset = 'COVID19_ERP127339', .keep = 'none') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/COVID19_ERP127339.csv')

deg_tibble <- read_csv('mission/SLE_TRPM2_MfMo/results/COVID19_ERP127339.csv')

deg_tibble |>
  mutate(adj.P.Val = padj, .keep = 'unused') |>
  write_csv('mission/SLE_TRPM2_MfMo/results/COVID19_ERP127339.csv')

dds_fpm <- fpm(dds)

dds_fpm |>
  as_tibble(rownames = 'gene') |>
  filter(gene == 'TRPM2') |>
  pivot_longer(-1, names_to = 'id1') |>
  left_join(sample_info) |>
  ggplot(aes(clinic, log1p(value), color = clinic)) +
  #geom_boxplot(outliers = F) +
  stat_mean(geom = 'col', fill = '#ffffff', alpha = 0) +
  stat_summary(geom = 'errorbar', fun.data = mean_cl_normal, width = .5) +
  geom_jitter(height = 0, width = .1) +
  labs(x = 'Disease severity', y = 'log(TPM+1)', color = 'Disease severity',
       title = 'TRPM2 expression in COVID19 patient PBMC',
       subtitle = 'PRJEB43380') +
  theme_bw() +
  scale_color_manual(values = c('royalblue','limegreen','darkorange','tomato'))
  
g1 <- last_plot()

write_csv(deg_tibble, '../covid19/results/wu-chen-deseq2-ALL.csv')

read_csv('../covid19/results/wu-chen-deseq2-ALL.csv') -> deg_tibble

deg_tibble %>%
  mutate(type = case_when(
    log2FoldChange > 1 & pvalue < 0.01 ~ 'Enriched in RR',
    log2FoldChange < -1 & pvalue < 0.01 ~ 'Depleted in RR',
    TRUE ~ 'NS'
  )) %>%
  mutate(sig_symbol = case_when(
    SYMBOL %in% deg_b_acti_list$SYMBOL ~ 'red',
    TRUE ~ 'grey'
  )) -> deg_tibble

ggplot(deg_b_acti_list, aes(log2FoldChange, -log10(pvalue))) +
  geom_point() +
  geom_vline(xintercept = c(1, -1), linetype = 'dashed') +
  geom_hline(yintercept = 2, linetype = 'dashed') +
  #scale_color_manual(values = c('red', 'blue', 'grey')) +
  ggpubr::theme_pubr()

# clusterProfiler -----------
deg_tibble %>%
  filter(log2FoldChange > 0, pvalue < 0.01) %>%
  pull(SYMBOL) -> up_list

deg_tibble %>%
  filter(log2FoldChange < 0, pvalue < 0.01) %>%
  pull(SYMBOL) -> down_list

clusterProfiler::enrichGO(up_list,
         OrgDb = 'org.Hs.eg.db',
         keyType = 'SYMBOL',
         ont = 'BP',
         universe = deg_tibble$SYMBOL,
         readable = TRUE) -> ora_go

ora_go@result %>%
  filter(p.adjust < 0.1) %>%
  separate(GeneRatio, into = c('count','size')) %>%
  mutate(GeneRatio = Count / as.numeric(size)) ->
  ora_go_res

ora_go_res %>%
  ggplot(aes(Description, GeneRatio, color = p.adjust, size = Count))+
  geom_point() +
  coord_flip() +
  ggpubr::labs_pubr() +
  scale_color_viridis_c()

clusterProfiler::dotplot(ora_go)

cnetplot(ora_go)

enrichplot::pairwise_termsim(ora_go) %>%
  enrichplot::treeplot()

# extract useful gene sets ---------
ora_go@geneSets %>%
  names() %>%
  clusterProfiler::go2term() -> term_collect

read_tsv('../covid19/results/wu-chen-go-pathways-digest.txt') -> digest_pathway

digest_pathway$Term <- reorder(digest_pathway$Term, digest_pathway$order, decreasing = TRUE)

digest_pathway %>%
  rename(count = Up) %>%
  rename(adjust.P = padj.up) %>%
  ggplot(aes(Term, gene_ratio, size = count, color = adjust.P)) +
  geom_point() +
  ggpubr::labs_pubr() +
  coord_flip() +
  expand_limits(y = 0)+
  scale_color_gradient(low = 'red', high = 'black')
